function [ data ] = addBackgroundModel(data)
load_settings;
disp('Add bg model: Calculating background model');
fns = cell(size(data)); % number of cameras

for ix = 1:size(data,1)
    for jx = 1:size(data,2)
        if ~isempty(data(ix,jx).cam_id)
            fns{ix, jx} = sprintf('%s/%08d/%08d_%06d.jpg', image_root,  data(ix,jx).cam_id, data(ix,jx).day, data(ix,jx).hour);
        end
    end
end

for ix = 1:size(data,1)
    disp(sprintf('Calculcating bg model for camera: %d', data(ix,1).cam_id));
    im = imread(fns{ix,1}); % Assumes there's at least one image per camera, which is a pretty safe assumption to make, unless something went horribly wrong
    im = rgb2gray(im);
    sizes = cell(1);
    sizes{1} = size(im);
    imStack = zeros([size(im), size(fns{ix}, 2)]);
    imStack(:,:,1) = im;
    for jx = 1:size(fns, 2);
        if ~isempty(fns{ix,jx})
            im = imread(fns{ix,jx});
            im = rgb2gray(im);
            sizes{jx} = size(im);
            im = imresize(im, sizes{1});
            imStack(:,:,jx) = im;
        end
    end
    [medianGradX medianGradY] = imageStackMedianGradients(imStack);
    [bg_im, ~]=reconsEdge3(medianGradX, medianGradY);
    
    % Future: parallelize this...
    for jx = 1:size(data,2)
        if ~isempty(data(ix,jx).cam_id)
            if size(bg_im) ~= sizes{jx}
                bg_im = imresize(bg_im, sizes{jx});
            end
            % calc hog of bg for each rectangle, append
            new_hog = zeros(size(data(ix,jx).ppl_hogs));
            for kx = 1:size(data(ix,jx).ppl_rects,1)
                cur_rect = data(ix,jx).ppl_rects(kx,:);
                bgm_rect = bg_im(cur_rect(2):cur_rect(2)+cur_rect(4), cur_rect(1):cur_rect(1)+cur_rect(3));
                bgm_rect = imresize(bgm_rect, reshape_size);
                hh = vl_hog(single(bgm_rect), binsize);
                new_hog(kx, :) = hh(:);
            end
            data(ix,jx).ppl_hogs = [data(ix,jx).ppl_hogs new_hog];
            
            new_hog = zeros(size(data(ix,jx).bg_hogs));
            for kx = 1:size(data(ix,jx).bg_rects,1)
                cur_rect = data(ix,jx).bg_rects(kx,:);
                bgm_rect = bg_im(cur_rect(2):cur_rect(2)+cur_rect(4), cur_rect(1):cur_rect(1)+cur_rect(3));
                bgm_rect = imresize(bgm_rect, reshape_size);
                hh = vl_hog(single(bgm_rect), binsize);
                new_hog(kx, :) = hh(:);
            end
            data(ix,jx).bg_hogs = [data(ix,jx).bg_hogs new_hog];
        end
    end
end


%% For each rectangle in data, compute the HOG feature of the background
% and append it

disp('Add bg model: Done with backgound model');
end